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Abstract. We obtain new closed-form pricing formulas for contingent 
claims when the asset follows a Dupire-type local volatility model. To 
5-H obtain the formulas we use the Dyson- Taylor commutator method that 

we have recently developed in [5jj6 , 8 for short-time asymptotic expan- 
sions of heat kernels, and obtain a family of general closed-form ap- 
I proximate solutions for both the pricing kernel and derivative price. A 

bootstrap scheme allows us to extend our method to large time. We also 
perform analytic as well as a numerical error analysis, and compare our 
i i results to other known methods. 

CM 

d 

l£3 Contents 

O h 1. Introduction Q] 

2. Theoretical Framework [7] 

CN 3. Closed Form Approximate Solutions [T7] 

4. Comparison and Performance of the Method [18] 

q 5. Option pricing with long maturity: the bootstrap scheme [22] 

CO References [27] 
(N 

d 

Q\ 1. Introduction 

O 

Financial derivatives (also known as contingent claims) are now a ubiqui- 
tous tool in risk management with approximately 600 trillion dollars worth 
^ of such contracts currently in the market. The pricing of such derivatives is 

therefore an active area of research in both Mathematics and Finance (see 



for example 12, 15, 17, 21, 32] and the references therein). In this paper, 
we will apply the perturbative (asymptotic) method introduced in [8] for 
numerically solving parabolic equations and then use this method to price 
European options. 

One of the earliest models used in pricing derivatives is the Black-Scholes- 
Merton model [3 ,27, for which the movement in the price Xt of the under- 
lying asset on which the claim is based is modeled by geometric Brownian 
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motion. For the Black- Scholes-Merton as well as for other models given 
by stochastic differential equations, the pricing of European options can be 
reduced to the calculation of certain solutions of parabolic equations, ob- 
tained through Ito's Lemma (and the change of variables t <— T — t) in the 
backward Kolmogorov equation. The resulting equation is a Fokker-Planck 
equation, which is an equation of parabolic type. Fokker-Planck equations 
more generally have important applications in statistical mechanics and in 
probability (see for example the monographs [4j[l6j[30] ) . Given that the as- 
set price is always assumed positive, the Fokker-Planck equation is solved 
on the positive half-line. One difficulty in treating this type of the equation 
is that the coefficients of the Fokker-Planck operator typically vanish at the 
boundary, making the equation degenerate. 

For example, for the Black-Scholes-Merton model, the resulting Fokker- 
Planck equation is given by 



d t U(t, x) - LU(t, x) = 0, < t < T, x > 
U(0,x) = h(x), x>0, 



(1.1) 
where 

(1.2) L := \j 2 x 2 d 2 + rxd x — r, 

is the Black-Scholes operator, a degenerate elliptic operator, t is the time to 
expiry, and h is the so-called pay-off function. For a European Call option 
with strike K and expiry (or exercice) date T, the pay-off function h is given 
by the formula h(xx) = \xt — K\+ := maxjxr — K, 0}, where xt is the price 
of the underlying asset at time T. Above, a and r are constant parameters, 
representing respectively the volatility of the underlying asset, and the cur- 
rent interest rate. Since the operator is degenerate at the boundary x = 0, it 
can be shown that the solution automatically vanishes there and no explicit 
boundary condition need to be imposed. 

A popular model related to the Black-Scholes-Merton model is the CEV 
model [To]. In the CEV model, the operator L is the form 

(1.3) L(t, x) = L(x) = -a 2 x 2a d% + rxd x - r, 

where a, a, r are constant. Yet another popular model is Dupire's local 
volatility model, for which we allow the volatility to change with time: 

L(t, x) = L(x) = -a 2 (x, t)x 2 d 2 + rxd x — r. 

Except in special cases, such as the Black-Scholes-Merton equation above 
and when L has constant coefficients, very few exact solution formulas to the 



problem (1.1) are available. It is therefore important to devise fast, accu- 
rate approximate solution methods. The focus of this paper is on obtaining 
approximate solution methods that are fast and accurate by combining stan- 
dard numerical methods with the asymptotic techniques developped in [8]. 
Fast solution methods are crucial when calibrating unknown parameters, 



especially in the Baeysian inference framework. We hope to address this 
question in a forthcoming paper. 

In view of the above discussion, it is justified to study the forward initial- 



value problem (1.1) for the general case when L is an operator of the form: 
(1.4) L(t) := -a(t, x) 2 dl + b(t, x)d x + c(t, x). 

We therefore allow for variable coefficients in both space and time. We as- 
sume throughout that a{x) > 0, for x > and that the coefficients a, b, 
c are smooth functions. The perturbative method introduced in [8j for the 
study of parabolic equations in arbitrary dimensions was fully justified in 
the case when a, b, and c and all their derivatives are bounded, and are 
bounded away from zero: a(x) > 7 > 0. In this paper we complete the re- 
sults of [§] with explicit formulas for the ID case. Then we numerically test 
our formulas for the Black-Scholes-Merton and CEV models, obtaining an 
excellent agreement between our theoretical results and the numerical tests. 



Both the Black-Scholes-Merton model (1.2) and and the CEV model (1.3) 
are more general than the models considered in (§] in that their coefficients 
do not satisfy the assumptions of the paper, yet the numerical tests indi- 
cates that the results of that paper are still valid for the more general models 
considered here. This observation suggests that the theoretical framework 
of |8| is applicable in greater generality. We plan to study this point in a 
forthcoming paper. 

To explain our method, let us recall that, under certain conditions on 
the operator L and initial value h, described in details in the next section, 



there exists a smooth function Qt(x,y) such that the solution to (1.1) has 
the representation 



(1.5) U(t,x) = / g t (x,y)h(y)dy 



The kernel function Qt(x,y) in (1.5) is the fundamental solution or the so- 



called Green function for the problem (1.1). 



Remark 1.1. Given that Qt(x,y) arises in several different contexts, we will 
call the function Qt(x,y) the transition density kernel, pricing kernel, heat 
kernel, or Green function interchangeably, depending on the context in which 
the object arises. 

As mentioned above, except for some very special cases no explicit formu- 
las for Qt(x, y) or U(t, x) are available. For the Black-Scholes-Merton model, 
a change of variables reduces the PDE to a heat equation that can then be 
solved explicitly. Therefore, exact formulas for the kernel Q BSM and the 
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solution JJ exist, which we recall now for further reference: 



gB8M {Xt y) = «p(-rt) exp ( | Hx/y) + (r - a* /2)t 



(1.6) 



yV2^aH V 2<r 2 t J 

rOC 

U BSM (t, x) = / Gf SM (x, y)dy = xN(d-) + Ke~ rt N{d^ 
Jo 



where J\f(x) = f x —r=e z2 l 2 is the cumulative normal distribution function 
(cumulative Gaussian distribution function) and 

d ± = Hx/K) + (r±a 2 /2)t _ 

However, for the time-dependent Black-Scholes-Merton model, where a and 
r are time-dependent, or local volatility models in general, closed form solu- 
tions are generally given by series expansions and difficult to use in practice 
or are not known (see, for instance, [j"H[23j ). 

The method that we use in this paper is to give an approximate closed- 



form solution for the equation (1.1) by giving an approximate closed- form 
expansion for the Green's function Qt{x,y). Since our approximation of the 
Green's function is in terms of Gaussian-type integrals, it gives a closed-form 
for the approximate price of a European call option for any one-dimensional 



model where the operator L is given by (1.4). In fact, as an application, we 
give the prices and Greeks (that is, suitable derivatives) of a European call 
option and perform an error analysis in Section [4} 

There exists a vast literature on obtaining asymptotic expansions of the 
Green's function Q t (x,y) when t small and x is close to y, especially in 
the case that L is independent of time [§[20j[22j[25|[28}[35j[35j[36]. (See 
also (T}[l4}[T8}[26}[34j). Many of these methods are based on a geometric 
interpretation of the operator L (or at least its principal part) as a Laplace 
operator on curved space, and require computing the geodesies in this space, 
which very often must be done numerically. Other approaches are based on 
pseudo-differential calculus. In particular, Corielli, Foschi, and Pascucci [7] 
use a parametrix construction for the problem (flTlj) to obtain a closed-form 



approximate solution. We recently developed in [5j|6j|8] a complementary 
approach to computing short-time asymptotics for Qt, based on parabolic 
rescaling, Taylor's expansions of the coefficients, Duhamel's and Dyson's 
formulas, and exact commutator expansion. We called this method the 
Dyson-Taylor commutator method. Our method is more elementary and 
appears very stable in practical implementations. 

Let us fix a function z = z(x,y) with the properties that z(x,x) = x 
and all its derivatives are bounded. The function z will represent the base- 
point for a parabolic rescaling of the Green's function. Then our short-time 
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asymptotics give an expansion for the kernel in the form: 
G t (x, y) = gf ] (x, y; z) + t^+D/^M (x> y . z) 



(1- 



Gl n] (x,y;z) = Gf ] (x, y; z) + t^G? (x, y; z) + • • • + t n / 2 G^ (x,y;z) 



where Ql" is the sum on the first n terms of the expansion and represents 
the n-th order approximate kernel, while i( n+1 )/ 2 S^ 1 is the remainder. The 
first term, g|°' is given by a dilated Gaussian function 



(1.9) 



Gf\x,y) = Gl\x-y) 



exp(- 



\x - y\ 



^2tvra(0,z) 2 2ta(0,z) 2 '' 



The Dyson- Taylor commutator method yelds an explicit algorithm to com- 



pute the terms G t for any n, if L is an operator of the form (1.4) and 
corresponding analogs in higher dimension. 

More precisely, our main result in [8] is that for the local volatility operator 



(1.4), the n-th order approximate kernel has the form 



(1.10) 



r l ' 2 y k {z,z + 



y 



y 



), 



t l/2 ' £1/2 > t V £l/ 2 

where the functions (z, x, y) are algorithmically computable (recall that 
z = z(x, y)). In this paper we shall compute the functions for k = 0, 1, 2 
at an arbitrary basepoint z. The details, based on the Dyson- Taylor com- 
mutator method method, can be found in Section [2] and 2.1 We therefore 
obtain new closed form asymptotic expansions of the Green function for lo- 
cal volatility models. In particular, the first order asymptotic expansion at 
arbitrary z = z(x, y) is given by 



(1.11) gl 1] (x,y;z 



1 



^2irta(S,zY 



! + 3a(0,z)a'(0,z)-2b(0,z) ^ _ ^ 



a'(0,z) , ,q , v (x 
:{x-yf + {x-z)^- 



2a(0,z) 2 
y) 2 -ta(0,z) 2 



g 2ta(0,z) 2 



2ta(0,z) 3 ^ yj v ; ta(0,z) 3 

We provide an explicit formula for the second order expansion of the Green 
function at the end of Section [2] This algorithm can be implemented very 
efficiently at least in dimension 1 and for n small, n = 1, n = 2. The numer- 
ical tests in Section [4] show that already the second-order approximation is 
adequate for the Black-Scholes and CEV models. 

For each term G^ in the expansion of the Green function, let denote 
the corresponding term in the expansion of the solution, 



(1.12) 



L{W{t,x) 



Gf ] (x,y)h(y;K)dy. 



Then using (1.5) and (1.8), we arrive at the expansion of the value of the 



contingent claim, 



(1.13) 



U(t, x) = U [n] (t, x) + t ("+l)/2gN ( tj x ) 

f/N (t, x) = (t, x) + i 1 / 2 ^ (t, x) + tU^ (£, x) + • • • t*/ 2 «M (t, x) 
where 

(1.14) t( n+1 )/ 2 c[ n ](i,x) := t(™ +1 )/ 2 / Sf ] h(y]K) = U(t,x) - U^(t,x) 



o 



is the remainder term (or error) in the expansion of the solution. In 
[8] we have shown that the remainder can be controlled in exponentially 
weighted Sobolev norms, when the operator L is uniformly strongly elliptic. 
These bounds on the remainder imply that, in this case, the error made 
by replacing Qt with Q\ n ^ in (1.5) is of order i n / 2 globally in space, the 



expected optimal rate. In [5], we consider degenerate operators, the symbol 
of which is strongly elliptic with respect to some complete metric of bounded 
geometry. For example, the Black-Scholes and the SABR models fit into this 
framework. By contrast, the CEV model with < (3 < 1 does not fit into 
this framework. Our numerical tests indicate nevertheless that the error 
term has the same order in t even for the CEV model with (5 < 1. For 
pedagogical purposes and error analysis we will list all the details for the 
time-dependent Black-Scholes and CEV models, although our results are 
more general. 

In Section [4] we perform a numerical error analysis by computing both 
the numerical solution U and expansion and estimating the error 

(1.15) \U(t,x)-U [n] {t,x)\ 

pointwise for the basepoint z(x,y) = x, when n= 1,2. The error analysis 
is in good agreement with the theoretical results, even though the local 
volatility operators considered in this paper do not necessarily satisfy the 
assumptions on the coefficients of L needed to establish the analytic error 
estimates performed in 5,6,8|. 

In Section [4] we then perform an error analysis. For the Black-Scholes- 
Merton model, for which an exact solution formula is readily available, we 
compare the expansions at the basepoint z(x, y) = x with the exact solu- 
tion. (Note however, that numerical errors arise also in the calculation of 
exact solutions, due to round-off errors and other approximations.) For the 
CEV model, we compare the expansions with benchmark formulas in the 
literature, in particular the Hagan- Woodward implied volatility approxima- 
tion [19]. 

Given that the kernel approximation is asymptotic in time, it guarantees 
good error control a piori only for sufficiently small t. In Section [5j we shall 
introduce a bootstrap scheme to extend our method to arbitrary large time. 
This strategy is based on the evolutionary property of the solution operator 



to (1.1). By doing so, we show that the error is remarkedly reduced. As 
an application in portfolio management, we also compute the Greeks (or 
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hedging parameters) of a European call option and compare our approxi- 
mations with the true Black-Scholes Greeks in Section 14.21 and Section HJ 
These applications again underline the accuracy of our methods. 

Acknowledgements. The authors would like to thank Marco Avellaneda 
for valuable suggestions and comments on the manuscript, and Jim Gather al 
for useful discussions. Victor Nistor also gladly acknowledges support from 
the Max Planck Institute for Mathematics, where part of this work has been 
performed. 

2. Theoretical Framework 

We begin by recalling the Dyson- Taylor commutator method, which we 
introduced in [6j[8], to obtain small-time asymptotic expansions for the so- 
lution of the initial-value problem: 

d t U(t,x) - L(t)U(t,x) = 

^ 2 '^ U(0,x) = h(x). 



Throughout the paper, the operator L will be given by (1.4), and we will 
omit the explicit dependence of L and of its coefficients on x. In addition, 
we tacitly assume that all the coefficients of L are regular enough to carry 
our the manipulations described next. For a rigorous justification in the 
case L is not degenerate, we refer to (6j[8]- 



If there is a unique solution to the initial- value problem (2.1), then the 
linear operator that maps the initial data h to the solution U is well defined. 
We refer to such operator as the solution operator. For constant-coefficient 
second-order operators, Lq, the solution operator forms a semigroup, de- 
noted by e tL °, t > 0; that is, the solution operator has the following prop- 
erties: 

(i) e tLo \ = I. 

(ii) e tlLo e t2Lo = e (' 1+ * 2 ) L °, t 1} t 2 > 0. 

The same conclusion hold for variable-coefficient, but time-independet op- 



erators L, under some conditions, for instance if L is strongly elliptic 29 
(that is, a(x) > 7 > for all x). When L is a time-dependent operator, 
L = L(t), the solution operator is no more a semigroup, but under some ad- 
ditional mild conditions, forms an evolution system S(t\,t2) [6||24]. For an 
evolution system, property ^ is replaced by S(ti,t2)S(t2, £3) = S(ti,t^), 
if < t% < t 2 < t%. Following the notation set forth in the Introduction, we 
denote the kernel or Green's function of the solution operator to the problem 



Ml by Ql 



t ■ 

Our method relies heavily on the study of distribution kernels of the evo- 
lution operators defined by our Fokker-Planck operator, so a brief discussion 
of distribution kernels and of our conventions is in order. 



Remark 2.1. Given a linear operator T mapping smooth functions with 
compact support into distributions, there exist a distribution kernel kx 

7 



such that 
(2.2) 



Tu(x)= I k T (x,y)u(y)dy. 



The integral above is interpreted as the pairing between test functions and 
distributions. In this paper, we will be interested in the integral representa- 



tion ( 2.2 ) in the case that T is a smoothing operator, that is, an operator that 
maps compactly supported distributions into smooth functions. Then, the 
kernel is a smooth function, and the notation kx(x,y) is justified point- 
wise. (For a more detailed dicussion, see for example [33].) In this case, 
we will write T(x,y) to denote the kernel and in general, we shall 

identify an operator with its distribution kernel. Let / be a smooth function, 
then we denote the operators of multiplication by / also with /. Additon- 
ally, we notice that there is no confusion when writing fT or Tf since the 
distribution kernels of these operators are given by fT(x,y) = f(x)T(x,y) 
or Tf(x,y) = T(x,y)f(y). Similarly, there is no confusion when writing 
d x T(x, y), since the distribution kernel of d x T is {d x hr){x, y) = d x (kT(x, y)). 
However k T9x (x,y) = -d y k T (x,y). 

We now introduce parabolic rescaling, which is a basic tool used in this 
paper. Let z be a fixed, but arbitrary point in K and s > a parameter. 
Given a function f(t,x) we denote by 

(2.3) f s ' z (t,x) :=f(s 2 t,z + s(x-z)), 

the parabolic rescaling by s of the function f about (0, z). Thus h s,z (x) := 
h(z + s(x — z)) for a function that does not depend on t. We will refer to 
z as the basepoint for the rescaling. Similarly we define a rescaled operator 
L s ' z by 

(2.4) L s ' z (t, x) := la s ' z (t, xfd 2 x + sb s ' z (t, x)d x + s 2 c s ' z (t, x). 



If U solves the initial- value problem (2.1), then U s ' z solves the rescaled 
problem 

dtU'>*(t, x) - L s > z {t, x)U s > z {t, x) = 
^ U s ' z (0,x) = h s ' z (x) 

Consequently, the Green functions of the operator dt — L and of the rescaled 
operator dt — L s ' z are related by 

gf(x,y) = s- x g^ t {z + a -\x -z),z + s-^y - z)) 

= t (z + t *(x-z),Z + t 2(y - Z)), if S = *2. 

We now proceed to compute the Green's function G^"' z of the rescaled 



problem (2.5) when t = 1. In order to do so, we shall consider the Taylor 



expansion in s at s = of the rescaled operator L s ' z , given in equation (2.4), 



up to order n. By "Taylor expansion" we mean that we Taylor expand the 



coefficients of L s,z and group all terms of the same order in s. The operator 
L s ' z can then be written as follows 

n 

(2.7) L" = £j'L' k + s n+1 L%. 1 <jk,x), 

where L*' , ^(t, x) contains all the remainder terms from the Taylor expansion 
of the coefficients. 

In this paper, we concentrate on calculating explicitly the second-order 
approximation of the Green function of L. Hence, we fix n = 2 from now 
on. For notational convenience, we denote g'(t,x) = -g^g{t,x) and g(t,x) = 
T^g(t,x). Then the second-order Taylor expansion in s of f s,z at s = is 
given by 

(2.8) r z (t, x) = /(0, z) + s(x - z)f'(0, z) 

+ s 2 tf(0, z) + s 2 {x - z) 2 f"(0, z)/2 + s 3 r(s, t, x, z), 

with s 3 r(s, t, x, z) the remainder. Below a = a(0, z) and all the other func- 
tions are to be evaluated at (0, z), unless stated otherwise. We then readily 
have the second order Taylor expansion of L s ' z in s at s = 0: 

(2.9) LI := ^a 2 dl, L\ = L\{x) := ad \x - z)d 2 x + bd x , 
and, L| = L\ x + tL\ t , where 

(2.10) L z 2x := {a' 2 + aa"){x- z) 2 d 2 x /2 + b'{x- z)d x + c, L z 2 t :=aad 2 x . 
Hence 

L s ' z (t,x) = L Z + sL\{x) + s 2 {L z 2 , x (x) + tL z 2jt ) + s 3 L s 3 ' z (t,x), 
where L^ z {t,x) is the remainder term. 



Remark 2.2. Each L z k in (2.7) has polynomial coefficients of order k in (x — z) 
and of order < k/2 in t. In particular, Lq is a constant coefficient operator, 
for which the Green's function is computed explicitly in (2.28). Thus, in or- 
der for the expansion to capture the time dependence of the coefficients, the 
coefficient must be expanded at least to second order in s. Time-dependent 
corrections will therefore appear only at order s 2 = t in the expansion of 



Let Qi be the Green function of the parabolic problem (2.1), that is, the 
solution is given by U(t,x) = f Qt(x,y)h(y)dy =: (Q[ h)(t,x). 

We begin the approximation scheme for by decomposing L into a 
constant-coefficient, second-order operator Lq, for which we can explicitly 
compute the solution operator, and a remainder: 

(2.11) L(t) = L + V(t) 



where V(t) is a time-dependent, variable coefficient, second order operator. 
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By DuhameFs principle we then have 
(2.12) gL = e tLo + [ e^-^Vin^dn. 







Repeated applications of Duhamel's formula leads to a recursive represen- 
tation of G+ time-ordered expansion: 

(2.13) Qt = e tL ° + f e ( *- Tl)Lo y(ri)e TlLo dTi 
Jo 

+ ( r\^- T ^ Lo V{T 1 )e^- T2)Lo V{T 2 )e T2Lo df + ■ ■ ■ 
Jo Jo 



+ f r ... f Td 1 e(*- Tl ) Lo y(ri)e (ri - T2)Lo y(r 2 ) • • ■ V {T d )e TdL ° df 
Jo Jo Jo 
*t f r i r T d+i 

P (t-T 1 )L 0V ( \Jti-t 2 )L 



+ / / •••/ e^ L W{r 1 )e^-^ L W(r 2 )---V{T d+1 )g!; i df 
Jo Jo Jo 

where, for notational convenience, we have set dr^ ■ ■ ■ dr 2 dri = df. This 
expansion can be rigorously justified, at least in the case when L uni- 
formly strongly elliptic and all the coefficients of L and their derivatives 
are bounded. See [6j[8] for details. In the limit d — > oo, it yields an asymp- 
totic time-ordered series, also called a Dyson series, for the Green's function. 
The integer d stands for the iteration level in the time-ordered expansion, 
which at this point is distinct from the order n of the Taylor expansion of 
the operator L. For consistency we need n > d [8]. We set from now on 
d = n = 2. 

A similar formula holds for the Green's function g^ s ' z of the solution 



operator for the rescaled problem (2.5). We recall that it is enough to 



compute an approximate Green's function at t = 1 for the rescaled problem 



by (2.6). We now choose the operator Lq to be exactly the zeroth-order 



Taylor expansion of L s,z , given in ( |2.7[ ). Then: 

V s ' z (t) := L s ' z (t) -L = sL\{x) + s 2 L z 2 (t, x) + s 3 L s 3 ' z (t, x). 



and using (2.13) with d = n = 2 and t = 1 yields 

(2.14) gr = e L o + slf + s 2 (Z^ + Z z tX + 1Q + K'>*, 

where 

T{ = f e^~ T ^ L oLfe nL odT 1: 
Jo 

Z*! = f r e^-^oLfe^-^Lle^dTzdn, 
Jo Jo 

Jo 

Zf,t= f e^ L onL z 2T e^ L odn. 
Jo 



(2.15) 
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Even though we set t = 1, we still keep the t dependence explicit in X| t 
to emphasize this term comes from Taylor expansion in t. The term 1Z S ' Z 
in (2.14) contains all the higher order terms and will be included in the 
remainder. 

The approximation for the Green's function of the original problem 
(1.1) is now obtained as follows. Let 

(2.16) T s ' z (x, y) =: G (x, y- z) + s G x {x, y; z) + s 2 G 2 {x, y; z) 

be the distribution kernel of the operator e L o + slf + s 2 (lf 1 +X| X +X| t ). 
The desired second order approximation is then given by 

(2.17) gf\x, y) = t-V2 T ^*(z + (x - z)/V~t, z + (y - z)/V~t), 

where z = z(x,y) is an admissible function. In particular, the kernels G^ 
appearing in (1.8) are given by 

Gf 1 (x, y) := r l / 2 G n (z + (x - z)/V~t, z + (y - z)/Vt; z(x, y)). 

We thus need to compute the distribution kernels of the operators e o, If, 
Ifi, 3-% x , Z| t . In order to do so, we exploit the semigroup property of e tL ° 
to carry out explicitly the time integration in (2.15). Before we proceed, we 
introduce some useful notation. 

By [Ti,T2] := X1T2 — T2T1 = — [T2, Tl] we shall denote the commutator 
of two operators Ti and T2. Two operators Ti,T% commute if [Ti,T2] = 0. 
For any two operators Ti and T2, we define adr l (T2) by ad,T 1 {T2) ■= [Ti,T2], 
and, for any integer j, we define adj,^(T2) recursively by 

adtp (T2) := adT^adi^ 1 (T2)) ■ 

We next recall a Baker-Campbell-Hausdorff-type identity proved and used 
in this setting in [8] (note that the operators Ti are unbounded). Namely, 
for any 6 £ (0,1) and differential operator Q = Q(x,d) with polynomials 
coefficients in x, we have 

(2.18) e eL oQ = P ad (Q,9,x,z,d)e eL K 

where P a d(Q, 0, x, z, d) is a differential operator with polynomial coefficients 
in x given by 

(2.19) P ad (Q, 8, x, z, d) = Q + £ -a4g(Q). 

i=i l - 

In proving this formula, we use the fact that the series is actually a finite 
sum, as we show below. In particular, P can be explicitly computed. 
A simple calculation gives the following lemma. 

Lemma 2.3. Let L m be a second-order differential operator with polyno- 
mial coefficients of degree at most m. Then ad 3 L z(L m ) = for j > m. In 
particular, we have [Lq,L|J = 0, ad 2 L ^{L\) = 0, and ad^z^L?,^ = 0. 

Proof. The proof is a simple calculation. □ 
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We now proceed to compute the integrals in (2.15) 



f 1 e^- T1 ^ L oLfe TlL odTi = + (1 -ri)[Lg,LJ])e L odri 

Jo Jo 



(Lf + -m,Lf])e^, 



^1,1 



C r e( 1 - T ^ L °L z 1 e {Tl - T2)L °L z 1 e T2L °dT 2 dT 1 , 
Jo Jo 

f P{L\ + (1 - Tl )[L* ,m(Ll + (1 - r 2 )[L* , Ll\)e L »dT 2 d Tl 
Jo Jo 



X 2>x = C e^-^Ll^ny^dn 
Jo 



\Ll x + (1 - nW^LlJ + Q-pF-m, [L z ,Ll x ]])e L odr 
L% x + \[L z Q ,Ll jX ) + ±[L* ,[LlLl x ]])e L °o, 



L 2,t 



C e^ L o Tl L z 2 T e^ L od n = f 1 r x L\ T e L °d Tl = h 
Jo Jo 1 



2± e ■ 



Hence ( |2. 14 ) becomes 
(2.20) 
where 

(2.21) Q 2 



[l + sQ l + s 2 Q 2 )e L « +K S 



~^2,x + 2 [-^0) ^2,x\ + g [^0) [-^0) ^2,x\] + 2^2,*' 



and 7£. s ' z is again the error term as in (2.14). Therefore, we only need 



to compute the commutators in the above formula to get the second-order 
approximation of QV' ■ 

We recall that we agreed that all functions in the commutator formulas 
below are evaluated at (0,2;). Hence a = a(0,z), a' = a'(0,z), and so on. 
We have 



(2.22) m,Ll]=a(0,zy a '(0,z)d: 



,3„/ a 3 



Ll,LlY = a b a"dl 
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and hence 

Lt[L z , Lf] = a A a' 2 (x - z)d 5 x + ba 3 a'%, 

{ ' } [L z , L\\L\ = a 4 a' 2 {x - z)& x + (b + 3aa!) a 3 a'^. 

To compute the other commutators, we need the following lemma, which 
can be proved by induction using that [AB, C] = A[B,C] + [A,C]B. In 
particular, [d 2 , x — z] = 2d x and [d 2 , (x — z) 2 } = 2 + 4(x — z)d x . 

Lemma 2.4. For i,j> 1 integers we have 

min{ij} . . . . 

d x (x-zy U (x)= £ k\Q^ k )(x-zy- k dt k u(x)- 

We therefore have: 

£y = « 2 (a 2 + aa") (x - z)d 3 x , 
{ ' ' +a 2 (b' + a ,2 /2 + aa"/2)dl 

and hence 

(2.25) [L z ,[L z ,Ll x ]] =a\a' 2 + aa"),d x , 
so that finally 

(2.26) (Lf) 2 = (aa'(x - z)) 2 d A x + 2(a 2 a' 2 + aa'b){x - z)d 3 + (aa'b + b 2 )d 2 x . 

It follows that the approximation kernel of QY * is given by the applica- 
tions of a differential operator with polynomial coefficients to the Green's 
function of e* L ° . If is a smooth function, we denote by the convolution 
operator with (j), then C ( f,f(x) :=</>* f(x) = J 4>(x — y)f(y)dy, which shows 
that the distribution kernel of is C^,(x,y) = <p(x — y). It is immediate 
to check that 

(2.27) dxCtj, = Cs x <j), 

while Cfj,d x = —Cqa- By Remark 2.2 the distribution kernel of e L o is given 



by 

(2.28) e L Hx,y) = -L=zM- X ~f ), a = a(0,z), 

V2vra 2 2a^ 

and hence e L o is a convolution operator. 

Then, by flggrfr d k e L o(x, y) = H k (G)e L o(x, y), where 8 = ^ and 
are the (rescaled Hermite) polynomials satisfying Hq = 1 and fZfc + i(0) = 
— Oilfc(G) + H' k (Q)/a 2 . The polynomials iifc are easily computed by induc- 
tion as: 

1 3G 

Fi(e) = -e, # 2 (e) = e 2 -^, J f/ 3 (e) = -e 3 + ^, 

69 2 3 „ „ K 10 „ 3 156 



(2.29) tf 4( e) = e 4 -^ + 4, H 5 (e) = -e 5 + ^e 



tf 6 (e) = e 6 -^e 4 + ^e 2 -^ 

ar or ar 
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Using (2.29) we have 

G (x,y;z) = e L ° 



(2.30) 
and 



1 exp(- |x - y|2 ' 



V2vra 2 



2a 2 



(2.31) 



G^x, y; z) = {[L x + l - [L , (x, y) 

= (bd x + aa'{x - z)d 2 x + ^a 3 a'd x ^ e L °(x,y) 

bHt(e) + aa'{x - z)H 2 {Q) + ^a 3 a'H 3 (0) ) e L ° 



1 



v / 2ro 2 



3aa' - 26 



+ (x - z) 



(x — y) 2 — a 



2o 2 

2\ 1 



(x-y)- ^( x -y) 3 



We now carry out a similar calculation for the next (and last) term of our 
asymptotic expansion, namely 

(2.32) G 2 (x,y-z) = (Q 2 e L o)(x,y), 



with Q 2 given by Equation (2.21). We finally have 

G 2 (x, y; z) = Ql^ r + L% x + \\L%, L% x ] + \[Ll [Lg, L Z 2 J] 



1 



(2.33) +2 L i' 2 + i L f[ £ o,^] + ffilfM + §[^^f] 2 J 

where Pj are polynomials in x — z and x — y with coefficients given in terms 
of the values of the functions a, b, and c, and their derivatives, all evaluated 
at z = z(x,y), as follows 

P = c, P 1 = b'(x-z), 



1 



1, 



1, 



Pi 



a a + a b + a a /2 + b + a (x — z 



12, 



+ a (bo! + a + a" (a; — z) 2 ) 



(2.34) 



P3 = a(x — z)(a'b + + 2~ aa ' 2 ) 



,2 r 



Pi 



1 •> 11 3 j 3 /9 / \2 

-a a + 2a a |a-a() + -a — z) 



P 5 = ^a' 2 (x-z), P 6 = JaV 2 . 
z o 

In particular, we obtain the following explicit formula. 
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Example 2.5. For the CEV model given by Equation (1.3), we have a 
o~z a , a' = aaz a ~ 1 , b = rz a , z = z(x,y), and hence, 



G? EV (x,y;z 




e 2a 2 z 2c 



Let us introduce the time dependent Black- Scholes-Merton model to cor- 
respond to the operator 



(2.35) 



L 



1 



o~ (tj x + r(t)xd x - r(t), 



Thus the difference between the usual Black-Scholes-Merton model ( 1.2 ) and 



the time dependent Black-Scholes-Merton model (2.35) is that in the latter 



we allow a and r to depend on time. Then the asymptotic formula for the 
time dependent Black-Scholes-Merton model is obtained by setting a = 1 



in the Example 2.5 since that formula does not contain time derivatives of 
the coefficients. 

At this stage, we can allow the basepoint z to vary with x and y. In 
Section 2.1 below we compute the expansion for the basepoint z(x,y) = x 



and compare it in Section |4j Different choices of basepoints z may lead to 
more accurate and stable approximations. In future work, we plan to study 
how to optimize the choice of z. 

Definition 2.6. We call a function z = z(x,y) admissible if z(x,x) = x 
and all derivatives of z are bounded. 

In |6,8j, we rigorously prove error bounds for the remainder term in 13 ) 
in Sobolev spaces under the assumption that z be admissible (and all the 
coefficients of L, together with their derivatives, be bounded functions, and 
L be uniformly strongly elliptic). The function z = z(x, y) can be thus quite 
general. 

2.1. Kernel expansions at z = x. The choice z = x yields a simplified 
expression for the approximation, since certain terms disappear, and the 
approximation yields the price of a European call option in closed form. In 
fact, the convolution with the approximate Green's function can be evalu- 
ated exactly and the price of a European call option given in closed form. 
In particular there is no need for numerical quadrature in evaluating the 
integrals, thus improving the speed of our calculations. 



Example 2.7. By setting z = x in(2.31) and evaluating all coefficient func- 
tions at (0, x), we obtain the first-order correction to the rescaled kernel Q\ B Z 
in the form: 

„ . , x — y _(*-y) 2 ( 3aa' — 2b a' n9 \ 
(2.36) G 1 (x,y;z = x) = ^=M=e ( _( x _ y )2). 



7TGr 
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Example 2.8. Similarly, the second-order correction to the rescaled kernel 
Qi"' z is obtained in the form: 

1 (x-y) 2 f 1 



G 2 (x,y; z = x) 



V2na 2 



e"^^ • { ^a 6 a' 2 H 6 (G) 



(2.37) + %- U 2 a" + 4aa' 2 + 36a') H 4 (Q) 
6 



+ * (a 3 a" + 2a 2 b' + 2aa + 2aa'b + a' 2 a 2 + 2b 2 ) H 2 (@) + c 



where Q 



x—y 
a(0,x) 1 



, and Hq, H A , H 2 are given by (2.29). 



Example 2.9. For the time- dependent Black- Scholes-Merton equation, we 
have b(t, x) = rx, c(t, x) = —r and a(t, x) = a(t) x so that 

(2.38) 



3cr 2 — 2r a ( x — y x 



V 2VT(7 Z X Z 



2cr 2 \ ax 

where all coefficient functions are calculated at (0,x). 

Example 2.10. The second-order correction to the rescaled kernel is given 
by 



G* SM (x,y;z = x) 



1 

'2-rrax 



-e 2 



x 



<7 4 + 4<ro-(0) + 4r 2 + 4rcr 2 
8a^ 



(2 39) 

4<t(0)<t + 4r 2 - 16m 2 + 15a 4 2 12r - 29cr 2 4 cr 2 6 



+- 



-X z + 



24 



-x 4 + -r 



80- 2 

where X = (x — y ) / (ax) . 

Example 2.11. For the time- dependent CEV model, c(t,x) = —r and 
a(t, x) = a(t) x a , b(t, x) = rx with a(0) = a so that 



(2.40) 



G? EV (x,y;z = x) = *~ V e~^Y 



?>aa 2 x a ~ 1 - 2rx 1 ~ a 
2cr 



aax (x — y 
ax a 



and 

Example 2.12. For the second order correction we have 

(2.41) G^ EV {x, y;z = x) = {P 2 H 2 + P 4 H A + P 6 H 6 - r) —1 e~fe£ 

V2irax a 



where 



Ho 



x-y 



a 2 x 2a I a 2 x 2a ' 



Ha 



x — y \ G(x — yY 



a 2 x 2a 



+ 



a 6 x 6a a 4 X 4a 
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x-y 

a 2 x 2a 



15 fx — y\ 4 45 fx — y\ 2 15 



+ 



a 2 x 2a I a 2 x 2a / a 4 X 4a V (J 2 X 2a / CJ 6 X 6a 



and 



P 2 = - (a 4 a(2a - l)x 4a - 2 + 2a(ar + <r(0) + aar)i 2a + 2r 2 x 2 ) 



P 4 = -a 4 ax 4a (a 2 (a - l)x 2a ~ 2 + Aa 2 ax 2a ~ 2 + 3r) , P 6 = -^Vx 8 "" 2 
6 8 

Note that in the above two Examples for the CEV model, setting a = 1 
leads to the corresponding approximation for the BSM model. 

3. Closed Form Approximate Solutions 

In this section, we consider European call options. For European put 
options similar results can also be obtained, either directly from the defini- 
tion or by using put-call parity [32]. In what follows, we will work with the 
expansion obtained by setting z(x,y) = x as the basepoint. In this case, 
we are able to compute the integrals defining the approximate option price 
jjWi from in closed form, which bypasses the need for more computa- 
tionally intensive integration methods such as numerical quadrature, which 
are needed for more general basepoints z. 



by 



By (2.36) and (2.6), the first-order approximate Green's function in given 



(3.1) Gl 1] (x,y) 



/ 2Tria 



2a 2 1 



3aa' — 26 . 



2aH 



(x - yf 



Similarly, by (2.37) and (2.6) the second-order approximate Green's func- 
tion is given by 

We recall that we implicitly assume all coefficients are evaluated at (0, x). 



(3.2) G l2] (x,y) 



'2ira 



(P 6 H 6 (~) + P 4 #4(H) + P 2 H 2 (~) - r) 



where S = the functions Hq,H4,H2 are given by (2.29), and 



(3.3) 



-Tg — —a a, , -T4 



[a 2 a" + Aaa' 2 + 36a') 



p 2 = - ( a 3 a " + 2a 2 6' + 2aa + 2aa'6 + a'V + 26 2 ) 



All the coefficient functions are evaluated at t = 0, z = x. 

For European Call options with strike price K, by (1.5) the n th - order 
approximated option price is 



(3.4) 



UW(t,x)= / gP(x,y)(y-K) + dy, 
Jo 
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where we only take n = 1,2 here. We recall that t = T — t is the time 
to expiry T and t is real time. So to be more precise, the n th -order option 
pricing formula for European call options with expiry time T is U^(T— t, x). 

We have already observed that the general form of the approximate ker- 
nel, when product of polynomial functions against a rescaled 



Gaussian. Therefore, the integration in (3.4) above can be carried out in 
terms of error functions. Explicitly, 

U [1] (t,x) = -^e~ (J ^f (2a-a'{x-K)) 



and 



+1 • (erf (%=5 N ) + I .)(/,/ + .,■- /v|. 



2 V V V2ta 



(3.5) 



U®{t,x) = uW{t,x) + T \ [erf (^T) + l) (K 

1 (*-y) 2 ft 2 n „ . ,n , . n 

+ . e • -a 2 a" - air + a' 2 12) t 2 

+ (t a + o"(z - K) 2 /3) t+-{tb 2 - a' 2 (x - K) 2 /6) 

a 

tba'(x-K) 2 a' 2 (x-Kf 
a 2 ' 4a 3 

Note that in financial applications (i.e., in the risk free measure) b(t, x) = rx 
and c(t, x) = —r. 

Example 3.1. For the CEV model we have 

rrfll , x O-X^y/t ix-K) 2 

U£ EV (t,x) = J. e 2**t**« ((2 - a)x + aK) 

(3 g) 

1 / „ I x — K 



and when a = 1, it reduces to the first order approximation for the Black- 
Scholes-Merton model. 

4. Comparison and Performance of the Method 

In this section, we discuss the accuracy and efficiency of our approxi- 
mation for the Black-Scholes-Merton and the CEV model. We employ the 
Black-Scholes-Merton model primarily as a didactic example, given that an 
exact kernel and option pricing formulas exists. For the CEV model, we 
compare our approximation to other solution formulas considered a bench- 
mark in the literature, in particular the Hagan- Woodward scheme |19| . 

What we find in general is a very good agreement of the approximate 
pricing formulas we derive in this paper with those available in the literature, 
but with significant advantage in the computational efficiency. In particular, 
the agreement is good even for times that are not small. In Section [5j 
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Figure 1 . Comparison of our first order approximation with 
the Black-Scholes formula. Parameters: K = 15, a = 
0.3, r = 0.1. Basepoint z = x. The left picture is for t = 0.1, 
and the right one is for t = 0.8. Note that the x-axis is scaled 
by 10, i.e. the label 150 means the stock price is 15. 

we propose a bootstrap scheme in time to improve the accuracy of our 
approximation for large time. 

4.1. Performance of the method. We start by discussing the Black- 
Scholes-Merton model, and choose the parameters K = 15, a = 0.3 and 
r = 0.1, and plot the exact and approximate solutions for < x < 25. 
We compare our formula with the Black Scholes exact solution formula for 
different times t. Figure [I] gives two different cases, which show that when 
t is small the two solutions are in very good agreement with an absolute 
error of order O(10 -3 ). We notice that even when t is not small, the error 
is small. Tables [T] and [2] give a analysis of the pointwise error for the first 
order approximation with respect to the exact Black-Scholes formula. 

Remark 4.1. Throughtout this section, we fix the basepoint z = x, so that we 
have closed- form approximate solution formulas, and we can better gauge 
the error introduced by the our method. For more general basepoints z, 
further error is introduced by the numerical quadrature used for the inte- 
gration and the truncation of the pay-off function h at large x (this error 
is lower order, however, if h is truncated at x large enough with respect to 
K). 

Remark 4.2. Formula (??) shows that the first-order approximation of the 
kernel depends linearly on rt. Therefore, the error grows more rapidly for r 
large at comparable times. The same observation holds for the CEV model. 
For Black-Scholes, this issue does not arise, since a change of variables allows 
to reduce to the case r = in the equation. 

Analytic pricing formulas for the CEV model in terms of Bessel function 
series have been derived for any value of /3 [9,13]. However, sum such series to 
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X 

t 


12 


13 


14 


15 


16 


17 


18 


0.01 


0.0000 


0.0000 


0.0313 


0.3266 


0.0387 


0.0019 


0.0000 


0.05 


0.0461 


0.3385 


0.0179 


0.3915 


0.0179 


0.4068 


0.3957 


0.1 


0.7 


0.7 


0.2 


0.5 


0.2 


0.4 


1.2 


0.2 


2.2 


0.3 


0.7 


0.9 


0.7 


0.3 


1.3 


0.5 


1.2 


2.1 


2.5 


2.7 


2.7 


2.9 


1.9 



Table 1. Error of the first order approximation for the BSM 
model, K = 15, a = 0.3, r = 0, error scale= 10 -3 . 



t 


12 


13 


14 


15 


16 


17 


18 


0.01 


0.0000 


0.0000 


0.1000 


0.0000 


0.9000 


2.0000 


3.0000 


0.05 


0.1 


0.9 


1.4 


0.1 


3.6 


8.7 


14.5 


0.1 


1.7 


3.8 


3.3 


0.3 


7.0 


15.9 


26.4 


0.2 


9.3 


10.7 


7.1 


1.2 


14.0 


30.2 


48.8 


0.5 


39.0 


31.4 


15.1 


8.4 


39.4 


76.0 


116.8 



Table 2. Error of the first order approximation for the BSM 
model, K = 15, a = 0.3, r = 0.1, error scale=10 -3 . 



accurate order can be very computationally intensive (but see Schroder |31| 
for methods to compute the pricing formulas more efficiently). 

The numerical tests show our approach yields accurate pricing formulas 
that are, however, computationally much simpler. We choose /3 = \,K = 
15, a = 0.3, r = 0.1 for parameters. Schroder [31| derived the exact CEV 
solution when (3 = |. Figure [2] gives the comparison of our method and 
the true solution of the CEV model for this value of /3 for different times. 
Again, we plot the two solutions for < x < 25. 

Hagan and Woodward in [19] studied more general local volatility models, 
for which the stock price under the forward measure follows the SDE 

dF t = 7 (t)A( J F t )dW t) 

for some deterministic and suitably smooth functions 7 and A. CEV fits 
into this general model. 

Using a singular perturbation technique, Hagan and Woodward obtain a 
very accurate formula for the implied volatility for this model. In the CEV 
case, their implied volatility reads 

ob = -j^ I 1 + — (1 - P)(2 + /3) j + ^ /2(1 _fl J , 
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Figure 2. Comparison of our first order approximation 
with the exact formula for the CEV model derived in |31| . 
Parameters:/^ = = 15, a = 0.3, r = 0.1. Basepoint 

z = x. The first graph is plotted when t = 0.1, and the 
second is when t = 0.5. In the graphs the x-axis is scaled by 
10, i.e the lable 150 means the stock price is 15. 



where 



e 2r(i-/3)T_i e rT S + K 
a = a \j 2r(l-/3)T ' /= 2 • 

The approximate pricing formula is then obtained from the Black-Scholes 
formula by using erg as volatility. 

When = §, the CEV formula can be computed exactly |31|. In this 
case, Hagan and Woodward's approximation is shown by Corielli et al to 
be very accurate |7j. We therefore take this approximation as benchmark 
for comparison with our method. In the following numerical comparison, we 
choose p = 2/3, K = 20, r = 0.1, a = 0.3 and different times r = 0.3, 0.5 We 
compute the prices on the interval [0, 30], and divide it into 300 subintervals. 
Since the prices near the strike is of most interest for practitioners, we 
compare the methods near K = 20. Figure [3] gives the results, from which 
we see that our approximation is more accurate than the Hagan- Woodward 
approximation near the strike for different times. 

We remark that our method can in principle yield arbitrary accuracy in 



the small-time limit if more terms in the kernel expansion ( 1.8 ) are included. 
Furthermore, it allows to derive approximate solution formulas for even more 
general models than those of Hagan and Woodwards. 

4.2. The Greeks. In this part, we use the second-order approximate so- 
lution to compute the Greeks of a European call option. The Delta and 
Gamma of a call option, collectively known as the Greeks of the option, at 
the point x are calculated as 

u(t, x + dx) — u(t, x — dx) 
d6lta= 2dx ' 

21 




Figure 3. Comparison of our approximation with Hagan's 
results for the CEV model near strike. Parameters: (3 = 2/3, 
K = 20, r = 0.1, a = 0.3. Basepoint: z = x. The first graph 
is plotted when t = 0.3, and the second is when t = 0.5. Note 
that in the graphs the x-axis is scaled by 10, i.e. the label 
200 denotes that the stock price is 20. 



and 

u(t, x + dx) + u(t, x — dx) — 2u(t, x) 

gamma = — ; 

dx z 

respectively, where u(t, x) is the option price. Some methods, for example 
the Monte Carlo method, can price options accurately, but they are not 
efficient for obtaining good hedging parameters. We shall show that our ap- 
proximations not only give option prices, but also Greeks accurately. Again 
for didactic purposes, we choose the Black-Scholes-Merton model for which 
the Greeks can be computed exactly. 

Since we can price options in closed form (by choosing z = x), we can 
calculate the Greeks in closed form by simply differentiating the approximate 
pricing formula. However, again because of the complexity of these formulas, 
we will obtain the hedging parameters numerically. 

In the numerical experiment, we choose the parameters as follows: ma- 
turity r = 0.5, volatility a = 0.5, strike K = 20, interest rate r = 10% In 
Figure |4j we plot the difference between our approximation and the exact 
solution for Delta when the stock price varies from to 40. Figure [5] does the 
same for Gamma. The numerical test shows that the pointwise difference is 
very small, of the order of 10 -3 in both cases. More specifically, the biggest 
error is around 13 x 10 -3 . 



5. Option pricing with long maturity: the bootstrap scheme 

The Dyson- Taylor commutator method gives an asymptotic expansion of 
the Green function in the limit t — > 0. Therefore, its accuracy is in priciple 
limited to times to maturity t relatively small. For long maturity options, 
we expect the error to be possibly large. In this section, we shall introduce a 
bootstrap strategy to price options with a long maturity time. The scheme 
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Figure 4. Comparison of Delta of our approach and the true 
values under the Black-Scholes model. Model parameters: 
t = 0.5, K = 20,a = 0.5, r = 10%. Basepoint: z = x. The 
left graph plots the delta computed by our method and the 
true delta. The right graph plots their difference. Note that 
in the second figure the scale is 10~ 3 . The x-axis is scaled by 
10, that is, the label 400 means the stock price is 40. 




Figure 5. Comparison of Gamma of our approach and the 
true values under the Black-Scholes model. Model parame- 
ters: t = 0.5, K = 20, a = 0.5, r = 10%. Basepoint: z = x. 
The left graph plots the gamma computed by our method 
and the true gamma. The right graph plots their difference. 
Note that in the second figure the scale is 10 -3 . The x-axis 
is scaled by 10, that is, the label 400 means the stock price 
is 40. 



is based on the properties of the solution operator. Let us illustrate the 
bootstrap in the time independent case. In this case, we recall that the 
solution operator forms a semigroup. The semigroup property then gives 
that 

(5.1) e tL = (e" L ) n , Vra G N. 
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Then, if n is sufficiently large, e« will be accurately approximated by our 
method. 

We next describe the bootstrap scheme, which can be rigorously justified 
at least for the case of strongly elliptic operators (a bounded away from 

zero) by the error analysis in |8|. In the bootstrap scheme, we use (g^\ 

to approximate e tL , where as before we denote the approximate solution 
operator by its kernel Q\ . Suppose now Q\j n is the second order approxi- 
mation, then the error is in the order 0((^ 3 ^ 2 ). Because there are n steps 
in the bootstrap scheme, the total error is in the order of 

0((t/nf/ 2 ) x n = 0(t 3/2 /V^), 

and consequently, for t fixed, it becomes smaller and smaller as n increases. 
A similar analysis shows that the bootstrap strategy with the first order 
approximation does not improve accuracy, given that in this case the error 
at each step is 0(t/n), so the total error after n steps is 

0(t/n) x n = 0(t), 

which does not converge to zero as n — > oo. 

We numerically tested this scheme for both the option prices and the 
Greeks. In the bootstrap scheme, closed-form approximate solutions are 
not available after the first time step, since we integrate the aproximate 
Green's function against an expression of the form (3.5), which contains 
error functions. Therefore, we must integrate numerically and introduce 
an additional error due to the numerical quadrature. This error can be 
controlled and made lower-order by choosing the space discretization step 
small enough. A further error, which can also be made lower-order, comes 
from the truncation of the integration at large x. 

In the first simulation, we used the Black-Scholes-Merton model, and set 
the parameters as time to maturity t = 1 (one year), strike K = 20, risk- 
free interest rate r = 10%, and volatility a = 0.5. The left graph of Figure 
^ displays the error of the first-order-closed form solution, the second- 
order closed-form solution, the first-order approximation with bootstrap, 
and the second-order approximation with bootstrap. We truncate the half 
line (0, +oo) at 200, and fix the number of the bootstrap steps as n = 10, that 
is the time step is At = 0.1. We choose the space discretization Ax = 0.1. 
From the graphs, it is clear that the second-order approximation greatly 
improves the accuracy compared with the first-order approximation. The 
bootstrap scheme with the second-order approximation reduces the error 
even further as expected (See Table [3] for a quantitative error analysis). 
As predicted, on the other hand the bootstrap scheme with the first order 
approximation introduces an extra error. 

We also notice that around S = 40 (the label 400 in the graphs) the 
error with the second-order bootstrap tends to increase, an effect of the 
truncation error. To verify it, we truncate the half line at 400. The right 
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Figure 6. Comparison of our first and second order approx- 
imation with or without bootstrap. Each curve is the differ- 
ence with respect to the Black-Scholes formula. Parameters: 
t = 1, K = 20, r = 10%, a = 0.5. Basepoint: z = x. In our 
numerical integration, we truncate (0, oo) at 200 for the left 
picture, and at 400 for the right one. Note that the x-axis 
scaled by 10, i.e. the label 400 means the stock price is 40. 



t 


3 


2 


1 


0.5 


0.2 


0.1 


error 


0.0268 


0.0379 


0.0177 


0.0038 


4.3682e-004 


3.5703e-005 



Table 3. Errors of the bootstrap scheme for different times 
under the Black-Scholes-Merton model. Number of boot- 
strap steps is fixed as 10. Parameters: K = 20, a = 0.5, r = 
10%. Errors are measured in L°°(0,40), and the benchmark 
is the Black-Scholes formula 



graph of Figure ([5| shows that the error does not tend to increase. We also 
tested the cases when the time to maturity is two and five years, obtaining 
similar results. 

To give a sense on how accurate our bootstrap scheme with the second 
order approximation is for large time to maturity, we repeat previous numer- 
ical simulation for difference times, and measure the error in L°°(0, 40). We 
recall that the trike price is at S = 20, we are taking a symmetric interval 
around it, and that the number of bootstrap steps is fixed at 10. We report 
the errors in Table [HI 

As predicted, we can increase the number of bootstrap steps to obtain 
arbitrary accuracy in the aproximation. Furthermore, for relatively large 
t the number of bootstrap steps should be correspondingly large, so that 
the compound error from each bootstrap step is under control at the end. 
For example, in our tests when t > 4, it is not enough to reduce the error 
by bootstrapping with 10 steps. Using only 10 steps in this case, in fact, 
introduces additional errors. For more detail, see [6]. 
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Figure 8. comparison of our first order and 
Figure 7. Butterfly second order approximation with or without 

option payoff. bootstrap for a butterfly option when t = 

l,r = 10%, a = 0.5. In our numerical inte- 
gration we truncate the half line at 200. Note 
that the x-axis is scaled by 10, i.e. the label 
400 means the stock price is 40. 



In order to eliminate the effect of the truncation error, we shall work 
with a butterfly option in the rest of this section. Mathematically, a butter- 
fly option corresponds to an intial pay-off given by a hat function, Figure 
0. Our method gives closed- form solution for butterfly options as well, 
by linearity. Figure ^ shows the errors of a butterfly option within the 
Black-Scholes-Merton model with K = 20, K\ = 15, = 25 obtained by 
our first order and second order approximation with or without bootstrap. 
The benchmark is the true solution. The parameters we were using are the 
same as we mentioned before. Again, we truncate the half line at 200. For 
a butterfly option, the truncation error is clearly very small, given that the 
data is compactly supported (Figure ([8])). For the second order approxima- 
tion with a bootstrap scheme, the error is almost zero. It is in the scale 
of 10~ 3 , while without the bootstrap the error is of the scale 10 -2 . This 
coincides with our theoretical results. 

We can also run the simulation as in Table [3j and we find comparable 
results. 

We conclude by discussing the bootstrap scheme for the Greeks. Directly 
using the closed-form approximation formula to compute the Greeks for 
very long maturity time (t » 1) is not advisable. In fact, our closed- 
form approximation for call option oscillates near the strike price, and the 
oscillation grows with the time to expiry, as the overall error grows. The 
appearance of the oscillation is due to the discontinuity of first derivative of 
the pay-off function at the strike price. This phenomenon is clearly visible for 
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Figure 9. Comparison of the delta computed by our boot- 
strap scheme ( 2nd order approximation) with the true val- 
ues under the Black-Scholes model. Model parameters: 
t = 0.5, K = 20,a = 0.5, r = 10%. Basepoint: z = x. The 
left graph plots the delta computed by our bootstrap method 
and the true delta. The right graph plots their difference. In 
our numerical integration, we truncate the half line at 400. 
Note that in the second figure the scale is 10 -4 . The x-axis 
scaled by 10, i.e. the label 400 means the stock price is 40. 



butterfly options, where the first derivative of the payoff function has three 
discontinuities, Figure ([8]). This oscillation is amplified in the calculation of 
Greeks. The bootstrap scheme reduces this oscillation dramatically. 

For the numerical simulation, we choose the same parameters as those in 
Section (4.2). The small time step ensures very good error control at each 
time step. Also, we minimize the truncation error as before by truncating 
the integral at x = 400 and comparing the approximations only on the 
interval [0, 40] near the strike price k = 20. The left graph of Figure ^ 
plots the true delta and our approximation in the same picture, and the 
right one plots the difference between these two curves, which shows that 
the difference between the true value and our approximation is in the order 
of 10 -4 with the biggest error around 3.3 x 10 -4 . Thus our approximation is 



quite accurate. For the gamma, we obtain similar results, see Figure 10 The 
difference is in the order of 10 -5 , and the biggest error is around 5.6 x 10~ 5 . 
In both cases, there are no oscillations on the same scale of the solution. 



References 

[1] Y. A" it Sahalia. Closed-form likelihood expansions for multivariate diffusions. Ann. 

Statist., 36(2):906-937, 2008. 
[2] R. Azencott. Asymptotic small time expansions for densities of diffusion processes. 

Lecture Notes Maths, 1059:402-498, 1984. 
[3] F Black and M. Scholes. The pricing of options and corporate liabilities. Journal of 

Political Economy, 81(3):637-54, May- June 1973. 

27 




Figure 10. Comparison of the gamma computed by our 
bootstrap scheme ( 2nd order approximation) with the true 
values under the Black-Scholes model. Model parameters: 
t = 0.5, K = 20, a = 0.5, r = 10%. Basepoint: z = x. 
The left graph plots the gamma computed by our bootstrap 
method and the true gamma. The right graph plots their 
difference. In our numerical integration, we truncate the half 
line at 400. Note that in the second figure the scale is 10~ 5 . 
The x-axis is scaled by 10, i.e. the label 400 means the stock 
price is 40. 



[4] H. Carmichael. Statistical methods in quantum optics. 1. Texts and Monographs in 
Physics. Springer- Verlag, Berlin, 1999. Master equations and Fokker-Planck equa- 
tions. 

[5] W. Cheng, R. Costantinescu, N. Costanzino, A. L. Mazzucato, and V. Nistor. Ap- 
proximate solutions to second order parabolic equations iii: the degenerate case. In 
preparation. 

[6] W. Cheng, A. L. Mazzucato, and V. Nistor. Approximate solutions to second order 
parabolic equations ii: time dependent case. Work in progress. 

[7] F. Corielli, P. Foschi, and A. Pascucci. Parametrix approximation of diffusion transi- 
tion densities. Preprint, 2009. 

[8] R. Costantinescu, N. Costanzino, A. L. Mazzucato, and V. Nistor. Approximate 
solutions to second order parabolic equations i: analytical estimates. Arxiv Preprint 
0910.1562v2, IMA Preprint 2248. Submitted. 

[9] J. Cox. Notes on option pricing 1, constant elasticity of diffusions, unpublished draft, 
Stanford University, 1975. 
[10] J. Cox and S. Ross. The valuation of options for alternative stochastic processes. 

Journal of Financial Economics, pages 145-166, 1976. 
[11] G. Dorfleitner, P. Schneider, K. Hawlitschek, and A. Buch. Pricing options with 
green's functions when volatility, interest rate and barriers depend on time. Quanti- 
tative Finance, 8(2):119-133, 2008. 
[12] D. Duffle. Dynamic asset pricing theory. University Press, 2001. 

[13] D. Emanuel and J. MacBeth. Further results on the constant elasticity of variance call 
option pricing model, the Journal of Financial and Quantitative Analysis, 17(4):533- 
554, 1982. 

[14] W. Farkas, N. Reich, and C. Schwab. Anisotropic stable Levy copula processes — 
analytical and numerical aspects. Math. Models Methods Appl. Set., 17(9): 1405-1443, 
2007. 

28 



[15] J.-P. Fouque, G.e Papanicolaou, and R. Sircar. Derivatives in financial markets with 
stochastic volatility. Cambridge University Press, Cambridge, 2000. 

[16] C. Gardiner. Handbook of stochastic methods for physics, chemistry and the natural 
sciences, volume 13 of Springer Series in Synergetics. Springer- Verlag, Berlin, third 
edition, 2004. 

[17] J. Gatheral. The volatility surface: a practitioner's guide. John Wiley and Sons, 2006. 

[18] P. Greiner. An asymptotic expansion for the heat equation. Arch. Rational Mech. 
Anal, 41:163-218, 1971. 

[19] P. Hagan and D. Woodward. Equivalent black volatilities. Applied Mathematical Fi- 
nance, 6(3):147 - 157, 1999. 

[20] E. Hsu. Stochastic analysis on manifolds, volume 38 of Graduate Studies in Mathe- 
matics. American Mathematical Society, Providence, RI, 2002. 

[21] J. Hull. Options, Futures and Other Derivatives. Prentice Hall, 2007. Sixth edition. 

[22] J. Kampen. The wkb-expansion of the fundamental solution of linear parabolic equa- 
tions and its applications. Submitted. 

[23] C. Lo, P. Yuen, and C. Hui. Constant elasticity of variance option pricing model with 
time-dependent parameters. Int. J. Theor. Appl. Finance, 3(4):661-674, 2000. 

[24] A Lunardi. Analytic semigroups and optimal regularity in parabolic problems. Progress 
in Nonlinear Differential Equations and their Applications, 16. Birkhauser Verlag, 
Basel, 1995. 

[25] H. McKean, Jr. and I. Singer. Curvature and the eigenvalues of the Laplacian. J. 

Differential Geometry, l(l):43-69, 1967. 
[26] R. Melrose. The Atiyah-Patodi-Singer index theorem, volume 4 of Research Notes in 

Mathematics. A K Peters Ltd., Wellesley, MA, 1993. 
[27] R. Merton. Theory of rational option pricing. Bell Journal of Economics, 4(1):141- 

183, Spring 1973. 

[28] S. Minakshisundaram and A. Pleijel. Some properties of the eigenfunctions of the 
Laplace-operator on Riemannian manifolds. Canadian J. Math., 1:242-256, 1949. 

[29] A. Pazy. Semigroups of linear operators and applications to partial differential equa- 
tions, volume 44 of Applied Mathematical Sciences. Springer- Verlag, New York, 1983. 

[30] H. Risken. The Fokker-Planck equation, volume 18 of Springer Series in Synergetics. 
Springer- Verlag, Berlin, second edition, 1989. Methods of solution and applications. 

[31] M. Schroder. Computing the constant elasticity of variance option pricing formula. 
the Journal of Finance, 44(1):211-219, 1989. 

[32] S. Shreve. Stochastic calculus for finance. II. Springer Finance. Springer- Verlag, New 
York, 2004. Continuous-time models. 

[33] M. Taylor. Pseudodifferential operators, volume 34 of Princeton Mathematical Series. 
Princeton University Press, Princeton, N.J., 1981. 

[34] M. Taylor. Partial differential equations. I, volume 115 of Applied Mathematical Sci- 
ences. Springer- Verlag, New York, 1996. Basic theory. 

[35] S. Varadhan. Diffusion processes in a small time interval. Comm. Pure Appl. Math., 
20:659-685, 1967. 

[36] D. Vassilevich. Heat kernel expansion: user's manual. Phys. Rep., 388(5-6) :279-360, 
2003. 

E-mail address: cheng@math.psu.edu 



Department of Mathematics, Pennsylvania State University, University Park, 
PA 16802 



E-mail address: costanzi@math.psu.edu 

29 



Department of Mathematics, Pennsylvania State University, University Park, 
PA 16802 

E-mail address: jcll2@psu.edu 

Department of Marketing, Smeal College of Business, and Department of 
Statistics, Pennsylvania State University, University Park, PA 16802 

E-mail address: mazzucat@math.psu.edu 

Department of Mathematics, Pennsylvania State University, University Park, 
PA 16802 

E-mail address: nistor@math.psu.edu 

Department of Mathematics, Pennsylvania State University, University Park, 
PA 16802 



30 



